Dynamic look-ahead feedrate scheduling method based on sliding mode velocity control

In the feedrate scheduling of complex curve direct interpolation, dynamic constraints such as axis acceleration and jerk are related to the actual state of the tool. Most existing methods convert dynamic constraints to velocity constraints at sampling points. However, it cannot guarantee the dynamic constraints are satisfied between sampling points. Addressing the issue, this paper proposes a dynamic look-ahead feedrate scheduling method based on sliding mode velocity control, which generates the motion command considering dynamic constraints in every interpolation cycle. To dynamically generate commands based on the current tool state, the acceleration and deceleration method based on sliding mode velocity control has been proposed, which can control tool state to transition to the command state with any initial state. To ensure sufficient distance for acceleration and deceleration, this paper uses braking distance to dynamically estimate the look-ahead distance. Then the minimum value within the look-ahead interval is selected as the command velocity for this scheduling cycle and the actual motion command is determined based on the dynamic constraints of each axis. Simulation and experiment results prove that compared with the existing method, this method effectively reduces the overshoot of dynamic constraints without significantly increasing the machining time. The analysis of real-time computation time has demonstrated the potential of the method proposed in this paper for real-time applications.


Principle of proposed method
The dynamic look-ahead feedrate scheduling method is shown in Fig. 1.The feedrate scheduling method is divided into two stages: non-real-time pre-processing and real-time tasks.The non-real-time pre-processing stage realizes the sampling of NURBS to obtain the differential properties of NURBS at sampling points.Combined with the feedrate constraints model such as chord error and axial velocity, the feedrate limitation profile satisfying the constraints is obtained.Meanwhile, the relationship between feedrate and look-ahead distance is established.The above information will be saved in the buffer to support the real-time tasks.
In the real-time task stage, the look-ahead distance of current period is calculated firstly according to the actual feedrate and the look-ahead distance function.After seeking the look-ahead interval, the minimum value of feedrate limitation profile in this interval is to serve as the command feedrate v cmd (k) of this period.Then SMVC method calculates the desired tangential acceleration and jerk command a cmd (k), J cmd (k) using the current tool status.Subsequently, dynamic constraints are used to adjust the a cmd (k) and J cmd (k) and the actual jerk command J act (k) satisfying dynamic constraints is derived.If there is no solution that meets the dynamic constraints, then J act (k) only needs to meet the maximum tangential jerk constraint.Once J act (k) is obtained, the arc length increment �s(k) and parameter increment �u(k) can be calculated.According to the �u(k) , the command position and NURBS differential properties of next period can be obtained.

Sliding mode velocity control acceleration and deceleration method
Define tool tip state x = v act a act T , where v act represents the tool tip's velocity and a act represents the tool tip's acceleration.The command state can be represented as x cmd = v cmd 0 T ,where the v cmd represents the com- mand tool tip's velocity.The A/D process is the process that tool tip's state transitions to command state x cmd from current state x .And in this process, the velocity and acceleration should be continuous and a act ∈ [−A max , A max ] .
Consider the following relationship between a cmd and v act where K 1 is a positive coefficient.The graph of Eq. ( 1) on phase plane is shown in Fig. 2. If a cmd satisfies Eq. ( 1), it will be continuous and in [−A max , A max ] .To prove the stability of control method described by Eq. ( 1), let v e = v cmd − v act and a e be the first-order derivative of the v e with respect to time t , which can shift command status from v cmd 0 T to 0 0 T .Equation (1) becomes If [v e , a e ] T can converge to [0,0] T , it implies that [v act , a act ] T can converge to v cmd 0 T .Given the maximum acceleration as J max , the K 1 should satisfy When the tool tip state satisfies Eq. ( 1) and K 1 satisfies Eq. ( 8), the system will converge to the command state x cmd .For states not satisfying Eq. ( 1), a sliding mode control method can be designed to make them approach the ideal state first.

Define Lyapunov function (1)
Define sliding mode surface as When s c = 0, the current tool state satisfies Eq. (1).When s c = 0 , it implies that the tool state is not on the sliding mode surface.By controlling s c to tend towards 0, the tool state can be ensured to approach the sliding mode surface.Consider following equation: where J act is the jerk of tool tip and K 2 is a positive coefficient.For the system described by Eq. ( 10), the Lyapu- nov function V 2 = s 2 c /2 is chosen.It can easily be proven that V 2 is positive definite and V2 is negative definite.It implies s c can stably converge to 0.
Increasing K 2 can improve the system's performance, but if K 2 is too large, it may cause chatter on the sliding mode surface.In this paper, K 2 is taken as 1/T s , where T s is the interpolation period.From Eq. ( 10), it can be deduced that (6) In the actual process, the jerk of the tool tip needs to satisfy the condition |J cmd | ≤ J max .It means Eq. (11) may not be satisfied.However, as long as J cmd doesn't change sign, the system's convergence direction will not change.The limit of jerk only affects convergence speed and the system remains stable.
In the A/D method, it is important to achieve the command velocity without overshoot.Near v cmd 0 T , as v cmd − v act approaches 0, Eq. (1) and Eq. ( 11) can be approximated as Then the system described by Eq. ( 11) can be approximated as a linear system as shown in Fig. 3.The transfer function of the system shown in Fig. 3 is The system depicted in Eq. ( 13) is a typical second-order system.The natural frequency and damping ratio are Using Cauchy's inequality, there is equality is achieved when Therefore, the actual system's damping ratio is greater than 1, indicating that the system is overdamped and will not overshoot during convergence.It should be noted that the system described by Eq. ( 13) is timecontinuous.Actual computer numerical control is a discrete system, which may cause some overshoot due to the discrete step size.
Figure 4 illustrates the process of adjusting the system state to the command state using SMVC, providing three initial conditions: Case 1: Initial velocity = 0, initial acceleration = 0; Case 2: 0 < initial velocity < command velocity, initial acceleration < 0; Case 3: Initial velocity > command velocity, initial acceleration > 0. It can be observed that in all 3 cases, SMVC can effectively control system convergence to the command state.
Figure 5 shows the response of SMVC to override changes.It shows that SMVC could flexible response the changes of override.
When the system enters the sliding mode surface, the convergence rate is determined by Eq. ( 1) and related to the velocity error.With the velocity error approaching 0, the convergence rate becomes very small, which can reduce the system's efficiency.Although Eq. ( 8) limit the value of K 1 , for specific command velocity, the system does not pass through the jerk extrema on Eq. (1) during convergence.So, without causing overshoot, K 1 can be appropriately increased to improve the system's performance.
Due to the nonlinear characteristics of SMVC, it is difficult to obtain an analytical expression between K 1 and overshoot.Therefore, this paper uses an iterative method to solve for the optimal K 1 .Let f (x) be the mapping from the parameter K 1 to the system's overshoot.The overshoot is calculated by v max − v cmd , where v max is the maximum velocity of the system during the given time response process, and v cmd is the command velocity.This mapping can be implemented through simulation.The process of using an iterative method to find the optimal K 1 is shown in Fig. 6.
For example, with v cmd = 50 mm/s, A max = 1000 mm/s 2 , J max =10,000 mm/s 3 , according to Eq. ( 8), K 1 = 0.0308 s/m is obtained.Given tolerance of 0.01 mm/s, a more optimal parameter K * 1 = 0.0385 s/mm can be achieved.Figure 7 presents the response of SMVC under K 1 and K * 1 , while also comparing them to the response of the (11) S-shaped A/D.It can be observed that the performance of the system after parameter optimization has become very close to S-shaped A/D method.At the end of entire curve, the command state is 0 0 T .If tool system directly follows the 0 0 T , it may cause velocity to be 0 before reaching endpoint.To address this issue, this paper introduces a position-errorfeedback velocity term in the endpoint command state, which changes the endpoint command state from 0 0 where p end is the endpoint position, p act is the current position, and K 3 is a positive coefficient.The added term ensures that the velocity of tool will not drop to 0 before reaching endpoint.

Feedrate scheduling method NURBS description
The definition of an p-degree NURBS curve is as follows: where N i,p (u), i = 0, 1, ..., n is the basis function of the p-degree B-spline defined on the node sequence , n are the control points of the NURBS curve, w i are the weight fac- tors corresponding to the control points P i , The B-spline basis function can be calculated by the de Boor-Cox recursive formulation:   For five-axis machine tools, machining curve is generally described by dual NURBS curve, which is composed of 2 NURBS curve, C(u) and (u) , as defined by Eq. ( 16).And C(u) represents the tool tip position and O(u) = H(u)−C(u) ||H(u)−C(u)|| represents the tool orientation vector.The motion relationship between five physical axes and the tool tip can be expressed as where q s (u) , q ss (u), q sss (u) are the first order, second order and third-order derivative of the axis position with respect to the arc length s .u s (u) , u ss (u), u sss (u) are the first order, second order and third-order derivative of the parameter u with respect to the arc length s.

Chord error constraint
To ensure the discrete points can reflect the characteristics of the original curve, it is necessary to consider the constraint of the chord error.
The chord error can be expressed as: where δ(u) is the chord error at the parameter u , ρ(u) is the radius of curvature of the curve at u , v t is the feedrate, T s is the interpolation period, and δ m is the maximum allowable chord error.The radius of curvature ρ(u) can be calculated by

Servo feed axis constraint
From Eq. ( 18), the velocity, acceleration, and jerk of the servo feed axis q can be expressed as: where v t , a t , J t are the tangential velocity, tangential acceleration and tangential jerk.Each servo feed axis has a maximum allowable velocity, acceleration and jerk, so it needs to satisfy:

Tracking error constraint
The tracking error of tool affect the machining quality and precision 11 , and limiting the tracking error of each axis could reduce the tracking error of tool.Currently, the servo feed axis usually adopts proportional-proportionalintegral control.And the control structure is shown in Fig. 8. where K p is the proportional gain coefficient of the position loop.K v is the proportional gain coefficient of the velocity loop.K vi is the integral gain coefficient of the velocity loop.M is the equivalent mass.B m is the equivalent damping coefficient and r g is the transmission coefficient.( 17)    q s (u) = q u (u)u s (u) q ss (u) = q uu u 2 s (u) + q u (u)u ss q sss (u) = q uuu (u)u 3 s (u) + 3q uu (u)u s (u)u ss (u) + q u (u)u sss (u) , q = X, Y , Z . . .The transfer function of the system is: Define the tracking error E(s) = X(s) − Y (s) , then: Equation ( 24) can be written in differential form: Where e(t), ė(t), ë(t), ... e (t) is the axis tracking error and its first-order, second-order and third-order derivative, J is the axis command jerk, a is the axis command acceleration, and v is the axis command velocity.According to the derivation of literature 26,27 , if Therefore, the tracking error of the axis is converted into a linear combination of velocity, acceleration, and jerk of the axis.The constraint of tracking error is: where

Feedrate limitation profile
The axis acceleration, jerk and tracking error are related to the actual system state.In this paper, the critical velocity satisfying above constraints can be estimated according to the state of a t = 0, J t = 0 , and the coefficient K dyn can be used to achieve a more conservative estimate.In this paper , K dyn = 0.9 .For satisfying the above constraints, there is: The feedrate limitation can be calculated as follows:

Dynamic lookahead distance
In the dynamic forward-looking feedrate scheduling based on SMVC, the look-ahead distance determines whether the system can satisfy constraints and smoothly follow command.To balance efficiency and smoothness, the look-ahead distance in this paper is calculated using the braking distance, which is the shortest distance required to decelerate to 0 at the current velocity and can obtained by numerical simulation of SMVC.Due to the constraints of acceleration and jerk, the actual braking distance may be longer than in the ideal case.Then a positive coefficient K 4 can be multiplied to ensure sufficient look-ahead distance.Since the look-ahead distance is velocity-dependent, command velocity oscillation may occur when approaching the velocity minimum, and jitter can be suppressed using a mean filter.

Interpolation algorithms
After obtaining the command jerk for this period, the state at the end of this interpolation period can be written: According to paper 28 , using the second-order Taylor development, there is:

Simulation
To prove the effectiveness of proposed method, simulations are carried out.The simulations were carried out on a PC with Windows10 operation system.The simulation environment was Matlab 2020a and the CPU was Intel® Xeon® Silver 4110 CPU @ 2.10GHz.The constraints and parameters for simulation of each axis are shown in Table 1:

Simulation 1: butterfly-shaped curve
A butterfly-shaped curve is used for simulation as shown in Fig. 9.And the curve information is shown in Appendix A Table A1.
The tangential feedrate command is set to 50 mm/s, the maximum tangential acceleration to 1000 mm/s 2 , maximum tangential jerk to 10,000 mm/s 3 , chord error limit to 0.1μm, K 1 = 0.0385 , K 3 = 5, K 4 = 2.The feedrate scheduling results in parameter domain are shown in Fig. 10a and results in time domain is shown in Fig. 10b-d, the axes tracking error is shown in Fig. 10e.The results show that the method proposed could schedule feedrate satisfying the tangential velocity, acceleration and jerk constraints.And the tracking error of each axis is also within the allowable range.
Axis acceleration and jerk are the typical dynamic constraints.Comparing the considering dynamic constraints time-optimal method 5 (M1) with the method proposed in this paper (M2), the results shows in Fig. 11 and Table 2.The results show that the method proposed in this paper, compared to M1, effectively reduces the exceeding of dynamic constraints without significantly increasing the machining time.

Simulation 2: open-pocket curve
The open-pocket curve is shown in Fig. 12.And the curve information is shown in Appendix A Table A2.
The simulation settings are same as Simulation1.And the feedrate scheduling results in parameter domain is shown in Fig. 13a and results in time domain is shown in Fig. 13b-d, the axis tracking error is shown in Fig. 13e,f.The tangential constraints and axes tracking error are all be satisfied.
Comparing M1 with M2, the results are shown in Fig. 14 and Table 3.Although M2 method has a greater number of exceeding constraints cycles compared to M1, the exceeding constraint ratio is very minor compared to M1.Therefore, it can still be considered that the method proposed in this paper has a better performance on satisfying dynamic constraints.

Experiment
To verify the real-time performance and tracking error constraints of the method proposed in this paper, experiments were carried out on a five-axis motion platform as shown in Fig. 15.
The control system of this platform is developed by TwinCAT3 and runs on a Windows 7 industrial computer an Inter® Core™ i7-6700 CPU @ 3.40 GHz.The interpolation and feedrate scheduling cycle is 2 ms.The tracking error coefficient can be identified through actual state and tracking error 4 .The identified result and constraints for each axis are shown in Table 4.

Experiment 1: butterfly-shaped curve with override change
The tangential feedrate command is set to 50 mm/s, maximum tangential acceleration to 1000 mm/s 2 , maximum tangential jerk to 10,000 mm/s 3 , chord error limit to 0.1μm, K 1 = 0.0385 , K 3 = 5, K 4 = 2. Before u ≥ 0.5, the over- ride is set to 100%.After u ≥ 0.5, the override is adjusted to 40%, which means the tangential feedrate command is reduced to 20 mm/s.After u ≥ 0.7, the override is set back to 100%.
The results are shown in Fig. 16. Figure 16a demonstrates the capability of the method proposed in this paper to adjust the override in real time.This capability provides the possibility for online control of cutting forces, which is an effective method for improving machining quality.Figure 16b shows that the maximum tracking errors for the X and Y axes are 0.4832mm and 0.4347 mm, both less than the constraint of 0.5 mm.This demonstrates the effectiveness of the tracking error constraints in this paper.and smaller tracking errors can indirectly reduce the machining contour errors, resulting in a surface with higher quality.

Experiment 2: open-pocket curve
Due to the dynamic performance of C axis, the experiment settings are same to experiment except K 4 =3.2.The results are shown in Fig. 17. Figure 17a demonstrates the capability of the method proposed in this paper for 5-axis machining.Figure 17b,c indicate that the maximum tracking errors for linear axes X, Y and Z are 0.0715mm, 0.4272 mm and 0.4636 mm, all less than the constraint 0.5 mm.The maximum tracking errors of rotary axes A and C are 0.0043 rad and 0.0114 rad, both less than the constrained value of 0.1 rad.

Computation time analysis
The computation time of Experiment1 and 2 is assessed by TwinCAT3 software and shown in Fig. 18.For the butterfly-shaped curve, the computation time can be kept below 50 μs in most cases, with a maximum not exceeding 80 μs.For the open-pocket curve, the computation time can be kept below 20 μs in most cases, with a maximum not exceeding 50 μs.Although it requires five-axis coordinate transformation and solution of properties of two NURBS curves, the computation time for the open-pocket curve is still significantly less than butterfly-shaped curve.This difference in computation time may be due to the differing complexities of the curves.The butterflyshaped curve has 52 control points, while the open-pocket curve has 8 (Supplementary information).
Figure 18a illustrates the variation in computation time under different override the override was reduced when u was within [0.5,0.7]. Figure 18a shows that after the reduction of override, the average computation time decreased.A possible reason is that after feedrate has decreased, the look-ahead distance becomes shorter, which means a shorter parameter look-ahead interval and fewer iterations required to find the minimum velocity.
Figure 18 demonstrates the computational efficiency of the method proposed in this paper for both complex planar tool paths and spatial five-axis tool paths, proving the great potential of the method for real-time applications.

Conclusion
High-precision direct interpolation of NURBS curves requires feedrate scheduling results satisfy the geometric, kinematic and dynamic constraints.Dynamic constraints are related to the current state of tool and cannot be directly constrained by velocity at sampling points.To achieve direct control of dynamic constraints, this paper proposes a dynamic look-ahead feedrate scheduling method based on sliding mode velocity control.The SMVC acceleration and deceleration method can control velocity to command velocity with any initial state and its stability has been proven.This paper also introduces the performance improvement method and end-pointreachable method.This paper analyzes the common constraints and linearizes the tracking error into a linear combination of velocity, acceleration, and jerk.The braking distance is used to estimate the look-ahead distance.The effectiveness of the proposed method in this paper is verified through simulation and experiment.In future research, selecting a smoother and more efficient sliding mode surface is a research direction The convergence speed at the endpoint also needs to be optimized, and further research is required on the selection of various coefficients.

Figure 1 .
Figure 1.Schematic diagram of dynamic look-ahead feedrate scheduling method based on SMVC.

Figure 10 .
Figure 10.Butterfly-shaped curve simulation results: (a) feedrate in parameter domain; (b) feedrate in time domain; (c) tangential acceleration in time domain; (d) tangential jerk in time domain; (e) tracking error.

Figure 13 .Figure 14 .Table 3 .Figure 15 .
Figure 13.Open-pocket curve simulation results: (a) feedrate in parameter domain; (b) feedrate in time domain; (c) tangential acceleration in time domain; (d) tangential jerk in time domain; (e) tracking error of X,Y and Z; (f) tracking error of A and C.

Figure 16 .
Figure 16.Experiment results of butter-shaped curve: (a) command and actual feedrate profile in time domain; (b) tracking error of axis.

Figure 17 .
Figure 17.Experiment results of open-pocket curve: (a) command and actual feedrate profile in time domain; (b) tracking error of X, Y and Z; (c) tracking error of A and C.

Figure 18 .
Figure 18.Computation time: (a) butterfly-shaped curve computation time and parameter timing diagram; (b) butterfly-shaped curve computation time histogram;(c) open-pocket curve computation time and parameter timing diagram; (d) open-pocket curve computation time histogram.

Table 1 .
Axis constraints and parameters for simulation.

Table 2 .
Performance on dynamic constraints of butterfly-shaped.